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A quantum-mechanical Gaussian wave-packet approach to 
the theoretical description of nuclear motions in a condensed- 
phase environment is developed. General expressions for the 
time-dependent reduced density matrix are given for a har- 
monic potential surface, and the exact quantum dynamics 
is found for a microscopic system-plus-bath model. Particu- 
lar attention is devoted to the influence of initial correlations 
between system and bath for the outcome of a pump-probe 
experiment. We show that the standard factorized prepa- 
ration, compared to a more realistic correlated preparation, 
leads to significantly different stimulated emission spectra at 
high temperatures. Recent experiments for the reaction cen- 
ter are analyzed using this formalism. 



I. INTRODUCTION 

The notion of wave packets is intimately connected 
with the foundations of quantum mechanics itself. De- 
spite of their importance during the initial stages of the 
development of quantum theory , this concept has been 
quickly overturned by the powerful and elegant opera- 
tor formalism. Only during the past decade interest in 
wave packets has emerged again, primarily triggered by 
the widespread availability of ultrafast femtosecond laser 
pulse techniques [||. Wave packets are by today a stan- 
dard tool employed to explain many different features in 
chemistry and physics, e.g., chemical reaction dynamics 
[^-^| , oscillatory motion of a coherent Bose-Einstein con- 
densate |6|, highly-excited Rydberg states in atoms 0, 
or electron-hole excitations in semiconductors ||||. The 
concept of wave packets is sometimes applied to ratio- 
nalize experimental data even though the reaction occurs 
under condensed-phase conditions, where the relevant re- 
action coordinate (e.g., describing the nuclear motion) of 
the wave packet ("system") may be strongly coupled to 
other modes ("environment" or "bath"). At this point 
one may ask whether it makes sense to use a wave-packet 
description for the reaction coordinate dynamics even if 
strong coupling to solvent modes is present. First, the 
description of the wave packet in terms of a wavefunc- 
tion is moot, and one has to use the reduced density 
matrix. Second, the bath leads to a damping of the wave 
packet and could cause the complete loss of coherence. 
In that case, the use of wave packets would be rather re- 
stricted. It is one of the purposes of this paper to clarify 
to what extent the wave-packet concept is applicable in 
the presence of strong system-bath coupling. 



The dissipation acting on the wave packet can have 
many different microscopic origins. In gas-phase reac- 
tions, one typically has rather weak damping due to cou- 
pling to the vacuum modes of the electromagnetic field 
(spontaneous emission), or due to collisions with other 
molecules. In contrast, dissipation can become decisive in 
condensed-phase reactions where one has strong coupling 
of the system to solvent or protein polarization modes. 
Considering previous theoretical treatments of gas-phase 
reactions, dissipation was mostly ignored or at best incor- 
porated within the framework of the Bloch equations |l(]] 
or the Redfield equations fill , where the latter allow to 
retain memory effects. However, as such an approach re- 
lies on perturbation theory in the system-bath coupling, 
its application to condensed-phase reactions character- 
ized by strong system-bath coupling remains question- 
able. Other methods are based on classical molecular 
dynamics (MD) simulations jl2| or projection operator 
techniques ||. 

In general, the problem of dissipative wave-packet mo- 
tion is theoretically quite demanding. In this paper, we 
treat a simple model introduced in Sections [n] and III 
but put particular emphasis on the effects due to dif- 
ferent initial states of the system arising in pump-probe 
spectroscopy experiments. This issue is shown to be im- 
portant for a correct description of experimental data 
on systems in a condensed-phase environment. More 
specifically, one might be tempted to assume a certain 
initial preparation which is named "factorized prepara- 
tion" henceforth. Under the factorized preparation, the 
density matrix at t = factorizes into a system part de- 
scribing the wave packet, and a part corresponding to 
the solvent modes. The wave packet at t = could then 
correspond to a pure state, e.g., a Gaussian wave packet. 
Here we provide a comparison of the factorized prepara- 
tion to the more realistic "correlated preparation" which 
takes into account initial system-bath correlations jb|. 
We stress that such preparation effects cannot be cap- 
tured by standard "dissipative wave-packet" approaches 
Pj , which implicitely use the factorized preparation. 

The structure of this paper is as follows. After present- 
ing general expressions for Gaussian wave-packet dynam- 
ics in Sec. |n|, the connection to microscopic system-plus- 



environment models is established in Sec. III. In Sec. I 



we then consider a pump-probe spectroscopy experiment 
involving two harmonic surfaces. As a practical example, 
we analyze the recent stimulated emission experiments of 
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the bacterial photo-synthetic reaction center by Vos et al. 
fbi] pHI , albeit the theory is more generally applicable. 
In that section, we also show that the two initial prepa- 
rations mentioned above cause pronounced differences in 
the emission spectra at high temperatures. Finally, some 
conclusions are offered in Sec. ^. Technical details have 
been deferred to an appendix. 



II. GAUSSIAN DENSITY MATRICES 

Let us start with general properties of the time evo- 
lution of a Gaussian reduced density matrix p(t). The 
Gaussian property implies that the underlying Hamilto- 
nian is at most quadratic in the system coordinate (?) 
and momentum (p), and imposes certain restrictions for 
the system-bath coupling. However, there is neither need 
to specify a Hamiltonian nor initial conditions for the 
system-bath complex at this stage [except consistency 
with the Gaussian form of pit)]. The spatial representa- 
tion of the density matrix reads 



p(q,q',t)=N- 1 (t)e^[-E(q,q',t)} , 



(2.1) 



where q',t) denotes a quadratic form, 
S(<?, q , t) = — g + —q -a 2 q-a 2 q+—qq (2.2) 

+ [a 2 + ^] 2 

a\+ a\ + 2a 3 

with arbitrary time-dependent coefficients ai(t). Further- 
more, the normalization Tr[p] = 1 is ensured by choosing 



N(t) = ^4Tr/[a 1 + al+2a 3 ] 



While a% = a' t + ia'{ and ai — a' 2 + ia' 2 ' can be complex- 
valued, p = implies a real- valued coefficient 03. There- 
fore we have five independent real-valued functions, and 
accordingly there are only five independent expectation 
values, 

(g) = 2a' 2 /[a' 1 +a 3 ] , 



(p) = h 



a-3 



([A?] 2 



0,3 
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([Aq,Ap] + ) = -h-^, 

where Aq = q — (q), Ap = p — (p), and [A, B} + = AB + 
BA. 

Let us now assume that the Ehrenfest theorem holds. 
This implies 



(p)/m 



-{[Aq} 2 ) = ([Aq,Ap} + )/r 



(2.3) 
(2.4) 



with the mass m. These equations eliminate two of the 
five degrees of freedom. Therefore we keep only (q(t)), 
the variance a(t) = 2([Aq] 2 (t)) , and the quantity 



Q(t) 



2h 



y/8a([Ap} 2 ) - m 2 a 2 



(2.5) 



as independent functions, where & = da/dt. The func- 
tions can be expressed in terms of these three quan- 
tities, and Eq. (|2 . 2f) then takes the form 



Re E ((g) 
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ImZ((q)+Q,(q)+Q',t) = 
&Q' d(q)\ (<jQ d(q) 



(2.6) 



(2.7) 



4rx 



dt J 



4rx 



dt 



The normalization constant becomes simply N = y/wa. 
The real-valued quantity is always within the bounds 
< 0(i) < 1, with the limiting case = 1 applying to a 
pure system. In fact, straightforward algebra yields 



6(i)=Tr[p 2 (i)] 



(2.8) 



It is noteworthy that 0(i) is in general an independent 
quantity, as there is no Ehrenfest relation expressing 
([Ap] 2 (<)) solely in terms of (q{t)) and a(t). 
By employing the unitary transformation 



U(t) = exp 



-mq 



11 
4cr 



d(q}_\ 
dt J 



exp 



-(q)p 



, (2.9) 



the density matrix p(t) 
independent form 



UpU> attains the coordinate- 



p(t) = Z- 1 exp(-/3i/) 



(2.10) 



with the Hamiltonian H(t) of a harmonic oscillator 
subject to an effective time-dependent confinement fre- 
quency 

h 



n(t) 



ma(t)Q(t) ' 
The effective inverse temperature f3 is 



(3{t) 



1 



In 
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(2.11) 



(2.12) 



and Z(t ) = [ 2sinh(fi»/3/2)1~ 1 . The transformed density 



Z(t) = \ 
rix ( PTO 



matrix ( [2.10 ) corresponds to the equilibrium density ma- 
trix of the harmonic oscillator H considered at fixed time 
t. 
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With the aid of this unitary transformation, it becomes 
easy to find the spectral decomposition of the density ma- 
trix. Transforming the result back to the original picture, 
we obtain 



P (t) = ^rx n (t)\<p n (t)){<p n (t)\ , 



where the eigenvalues are given by 
26 /1-6 



A»(f) = 



(2.13) 



(2.14) 



i + e vi + 

The spatial representation of the eigenfunctions is 

Vn(Q + (q),t) = ( 7 r ( re)- 1 / 4 (2"n!)- 1 / 2 J ff„(Q/N^e) (2.15) 
i m 3 &(d(q)/dt)' 2 



III. MICROSCOPIC SYSTEM-PLUS-BATH 
MODEL 

Let us now consider a wave packet moving in a har- 
monic potential surface under the influence of a bath 
composed of harmonic oscillators. The harmonic oscil- 
lator modes need not correspond to physical modes but 
could represent effective modes chosen to mimic the ac- 
tual environment in an optimal way. Such a system-plus- 
bath model allows us to derive exact expressions for the 
independent expectation values (q(t)), <j(t), and 6(i) of 
Sec. [n], an d thereby to obtain the exact quantum dy- 
namics of the damped wave packet for a given spectral 
density of the bath modes. In this section, the simpler 
case of a factorized preparation is treated. The correlated 



xexp,- 



x exp 



8([A P ] 2 ) 
Q 2 imQ 



2aG 



bQ_ d(q}_ 

4<7 dt 



where H n are the usual Hcrmitc polynomials. Trans- 
forming also H back into the original basis, we obtain 



H(t) = 



(Ap - m&Aq/2<j) 2 m0 2 [Aq] 



2m 



(2.16) 



In the end, the density operator can be written in the 
coordinate-independent form 



p{t) = Z- l exs>[-p(t)H(t) 



Using Eq. ( 2.14 ), one readily checks that 

oo 

T*Lo(*)] = X>„(*) = l, 



(2.17) 



preparation is then discussed in Sec. IV B 



We study a system-plus-bath model M, H = H$ + 
Hb + Hi, where the "system" part describing the un- 
damped coherent nuclear motion reads 



p 2 mto 2 , , 



(3.1) 



The "bath" is composed of harmonic oscillators coupled 
linearly to the system coordinate, 



Hr 



Hj 



9 ? 
Pj 1TljU}„ 
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(3.2) 

The influence of the bath onto the system is fully specified 
by the spectral density, 



j h = 77 y] — —s(u-uj) 



(3.3) 



and similarly one recovers Eq. (2.?), since 
Tr[p 2 (<)]=]TA 2 (i) = eW. 

n 

The linear entropy is then given by 

^ in (i) = l-Tr[p 2 ] = l-e(t), 

and the Shannon entropy is 

S(t) = -Tr[p\np] 
= A„ In A„ 



i, /i-e 2 \ 1 , fi + & 
2 ln lT + 2e ln lT^e 



(2.18) 



(2.19) 



A frequently used model spectral density is given by the 
ohmic bath with a Drude cutoff || , 



1 + {u/ljd) 2 



(3.4) 



Under the factorized preparation, the density matrix at 
time t — factorizes according to 



p(t = 0) = ps(q) ® Pb({xj}) 



(3.5) 



Here ps(q) describes a pure Gaussian wave packet of 
width (To centered around q = 0, corresponding to the 
wavefunction 

i/>(q) = (vr ( 7o)- 1/4 exp(-g 2 /2ao) . 

The bath is assumed to be in a thermal distribution, with 
the system coordinate held fixed at q = 0. While equi- 
librium properties of a damped harmonic oscillator have 
been studied exhaustively in the past, the consideration 
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of wave-packet initial preparations and their correspond- 
ing time evolution leaves room for our contribution. 

Due to the harmonic nature of the total system-plus- 
bath complex, the exact time-dependent density matrix 
p(q,q',t) can then be directly obtained from Feynman- 
Vernon theory Switching to symmetric and anti- 

symmetric linear combinations, 



x = (q + q')/2 , 



y 



(3.6) 



the propagating function of Ref. [[l3[ immediately leads 
to the result 

p(x, y,t) = -±= exp{~(x - (q)) 2 /a - y 2 /A<jQ 2 (3.7) 
\/7rcr l 



-my 



&(x - (q)) d(q) 
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dt 



}• 



in accordance with the general form (2.1). Now the three 



independent expectation values can be expressed in terms 
of microscopic parameters, 



(q(t)) = U$q / dt'G(t') , 
JO 

a(t) = a G 2 (t) + ^—G 2 (t) + ^K q (t) , 



(3.8) 
(3.9) 



and O(i) is defined by Eq. (2.5) with 



([Ap] 2 ) = — ( <J G 2 (t) 



m z <jQ m 



(3.10) 



The definition of the functions G(t), K q (t), and K p (t) 
for an arbitrary spectral density J(u>) can be found in 
the appendix. From these expressions, one verifies that 
the correct equilibrium values erg and ([Ap] 2 )^ M are 
approached at long times. 



Using the ohmic spectral density (3.4), we now discuss 
the question of coherence of the damped wave packet. 
Above a critical value 7 C of the damping strength 7, 
where j c follows from Eq. ( A14 ), oscillations in (q(t)) 
disappear and only incoherent relaxation can take place, 
see Figure |]. For u) D 3> u>o, the critical damping strength 
is given by -f c — 2ujq B. While this limit is of most in- 
terest in solid-state applications, the regime lod ~ ^0 as 
well as ujd <C loo has many applications in chemical sys- 
tems. Interestingly, the value of j c increases when u>£> /u>q 
becomes small. In fact, for u>d — > 0, the dynamics is al- 
ways fully coherent, 7 C = 00. In that limit, the bath is 
too slow to cause relaxational behavior. It is noteworthy 
that the precise value of 7 C depends on the quantity con- 
sidered in defining coherence. Taking the disappearance 
of the inelastic peaks in the spectral function as the rel- 
evant criterion leads to a critical value that is smaller by 
a factor l/y/2 [ p~7| . Since our coherence criterion is based 
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FIG. 1. Critical damping strength y c as a function of 
(solid curve). The limiting value y c /2uio ~ 1 for ujd 2> u)o is 
indicated by the dashed line. 

on the oscillatory behavior of (q(t)), where the latter 
equals the corresponding expression for a point-like par- 
ticle, the coherent-to-incoherent transition occurs at the 
same damping strength 7 C for a wave packet and a point- 
like particle. In particular, j c takes a temperature- 
independent value. 

Next we briefly discuss the time dependence of the vari- 
ance a(t), see Figure||], and of the Shannon entropy S(t), 
see Figure 0. The initial width ctq of the wave packet 




FIG. 2. Variance a(t) as a function of 7 for h{3uo — 1, 
o"o — 1-5 h/mujo, and ojb/ljo = 5. 

mainly influences the dynamics during the initial stage 
of the relaxation. Expanding for small times St, the vari- 
ance reads 



a (St) ~ <Tq + 



m 2 cro 



■ a (u 2 + ju>d) 



St 2 



(3.11) 



Therefore the variance initially increases (decreases) for 
(Jq < a (cr > cr)j where a = hjra^Juy^ + "/ujd- For 
7 < 7c oscillations in both a(t) and S(t) are found, 
similar to the behavior of (q(t)). These oscillations again 
persist at high temperatures, albeit with smaller ampli- 
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FIG. 3. Shannon entropy S(t) as a function of hficJo for 
7/2wo = 0.1, lud/loo — 15, and ao = 1- 

tude. The initial entropy increase observed in Fig. |^ be- 
comes very pronounced if ao strongly deviates from the 
natural width op of the damped oscillator. Furthermore, 
initial transient oscillations then persist for a longer time. 
They are particularly pronounced for low temperatures 
and small ao , with a transient entropy large compared to 
the equilibrium value S(t — > oo). 

IV. PUMP-PROBE SPECTROSCOPY OF THE 
REACTION CENTER 




q 

FIG. 4. Pump-probe setup involving two harmonic 
surfaces. The dark (excited) state surface is centered at 
Q = (q = go). 

make the (strictly speaking unphysical) assumption that 
the pump pulse transfers the complete nuclear wave 
packet up to the excited state surface. Therefore we only 
have to treat the dissipative excited state dynamics up to 
the probe pulse. Of course, thereby potentially important 
effects like the impulsive resonance Raman contribution 
P,pyT^] are missed. However, in principle our theory can 
straightforwardly be extended to provide a more realistic 
modelling of the pump and probe processes. 



Next we apply the results presented before in a specific 
context. The system under study is the photosynthetic 
reaction center in purple bacteria. In recent pump-probe 
experiments on modified and wild-type reaction centers, 
Vos et al. p^-p^| have observed oscillations in the time- 
resolved emission signal, which were interpreted to re- 
flect coherent nuclear motion in the excited electronic 
state ( "vibrational coherence" ) . This observation im- 
mediately received much attention, as coherent dynam- 
ics was not expected to exist in such a condensed-phase 
system. Clearly, a nuclear coordinate within a macro- 
molecule like the reaction center could be drastically in- 
fluenced by dissipation, which suggests a treatment sim- 
ilar to the one discussed above. The situation that Vos 
et al. constructed from their data is depicted in Fig. ^. 
The excited state surface was found to be parabolic with 
a curvature of ujq — 75 cm . It was populated with a 
870 nm pump pulse, say, at time t = 0, and probed with 
pulses around 921.5 nm, corresponding to the minimum 
of the excited state surface. In this section, we expand on 
the above analysis in order to describe the emission sig- 
nal. Thereby effects of the spectral density characteris- 
tics and of the initial correlations can be captured, where 
especially the latter are missed by any simpler formalism. 
In order to clearly show these initial correlation effects, 
we shall crudely simplify the modelling of the pump (and 
to a lesser extent of the probe) pulse. In particular, we 



A. Model and parameters 

The Hamiltonian H(t) = Ho+V(t) governing the emis- 
sion process first consists of an unperturbed Hamiltonian 

H Q = \G)H G (G\ + \E)(H E +H X + H B )(E\ . (4.1) 

The orthonormal states \G) and \E) denote the electronic 
degrees of freedom, with Hq (He) being the Hamiltonian 
in the ground (excited) state, 



H G = 



P 

2m 

P 



iG q 2 



2m 2 

where Tlloa = Eo — Eg and qo is the separation of the 
potential minima, see Fig. ||. The dissipation acting on 
the wave packet in the excited state is included via Hi + 
Hb, see Eq. (3.2). Since we only consider the dynamics 



on the excited state surface between the pump and the 
probe pulse, it is not necessary to account for dissipation 
in the ground state at t > 0. The effect of the probe pulse 
is described by V(t). Under the dipole and the rotating 
wave approximation [^J, 



V(t)=e(t)\G)(E\ + H.c. , 



(4.2) 
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where e(t) represents the temporal envelope of the electric 
field. The probe pulse was taken in the form 



;(£') = 0[t' - (t - S)} 0[t + 6- t'] 



(4.3) 



where 8 is the Heaviside function and the probe pulse 
is centered at time t. Since the 30 fs probe pulses used 
in Ref. |14| ] did not maintain their full intensity over the 
whole pulse duration, we have chosen a smaller duration 
of 28 = 20 fs. 

The model parameters were taken as follows. The 
frequency of the excited [ground] state surface is loq = 
75 cm -1 [log — 130 cm ]. Furthermore, qo and Hloa 
are calculated from the wavelength of the pump pulse, 
Ap« = 870 nm, and of the probe pulse acting at q = qo, 
Ap r = 921.5 nm. The parabolic geometry of Fig. || then 
yields qo = 2.07 'y^h/mcoo and loa = 11327 cm -1 . At 
this point, little is known about microscopic details of the 
dissipation acting on the reaction coordinate q. In prin- 
ciple, one should first compute the appropriate spectral 
density for the system under consideration by means of 
MD simulations (l!| . In the absence of such information, 
we make the assumption of an ohmic bath with a Drude 
cutoff, see Eq. (\j.4j). This spectral density was shown to 
be in agreement with the overall structure of the spectral 
density coupling to the primary electron transfer step in 
the reaction center 19 1 . To account for the lack of knowl- 
edge concerning the spectral density, we have studied two 
different spectral parameter sets. The first one, which is 
referred to as SP I, is j/2lo = 0.1 and lojj/loo = 100. 
The second one, referred to as SP II, is j/2lo = 0.75 
and lod/loq = 0.5. Both sets are within the coherent 
regime, 7 < 7 c (w.d), and are chosen such that the oscil- 
lations in (q(t)) decay on the same time scale as those of 
the T = 10 K emission signal reported in Ref. Jl4[ . 



B. Initial preparation 

A conceptually more severe point concerns the proper 
description of the initial state (t — 0). Again, two very 
different initial preparations are conceivable. The first 
one is to assume a wave packet in the usual sense, where 
the oscillator is initially in a pure state without correla- 
tions with the bath. This is the "factorized preparation" 
elaborated in Sec. Ill and (at least implicitely) employed 
in most previous treatments. On the other hand, the 
nuclear coordinate already experiences the environment 
while the system is in the ground state, and therefore the 
initial density matrix does not factorize. A more realis- 
tic preparation is to take the \G) oscillator at equilibrium 
with the same bath as in the excited state, whence there 
will be system-bath correlations at t — ("correlated 
preparation"). We mention in passing that the corre- 
lated preparation is related to the initial bath preparation 
discussed in Ref. [EGl in the context of electron transfer 



reactions. The special case log — w o with loo S> loq has 
also been treated in Ref. |^| and references therein. 

Due to its very short duration, as a result of the pump 
pulse at t = 0, the system is assumed to suddenly change 
from the ground state to the excited state surface. This 
amounts to both a vertical shift and a change in curvature 
log — * ^0, see Fig. [|. Technically speaking, the correlated 
preparation can be most conveniently accounted for by 
following the path- integral analysis of Ref. [jl3| , but keep- 
ing different system potentials acting on the imaginary- 
time and real-time paths. The resulting reduced density 
matrix is then of the form ( ft.Tp again. Due to the Ehren- 
fest theorem, (q(t)) and (p(t)) coincide with the results 
of the factorized preparation. The variances <j{t) and 
([Ap] 2 (<)) follow in closed form and are given in the ap- 
pendix. For the corresponding results under the factor- 
ized preparation, see Eqs. ( |3.9[ ) and ( 3.10 ). 

For both preparations, the initial width a(t — 0) was 
chosen as the thermal width cr? in the ground state os- 
cillator. Importantly, despite of having the same initial 
value, the time-dependence of the variance is strikingly 
different depending on the initial condition. This be- 
comes particularly evident for log = wo, where for the 
correlated preparation, a(t) and ([Ap] 2 (t)} stay constant 
in time, whereas the factorized preparation always leads 
to time-dependent variances. This can be understood 
by noting that ([Ap] 2 (i = 0)) for the factorized prepa- 
ration is determined by the minimum uncertainty con- 
dition <d(t = 0) = 1, see Eq. (|2.5|), while it is given by 



([Ap} 2 }9, see Eq. (A8), in the case of a correlated prepa- 



ration. Since the deviation in ([Ap] 2 (i = 0)) for the two 
initial preparations becomes larger with increasing tem- 
peratures, one expects that the choice of the correct ini- 
tial preparation is more important at high temperatures. 
This is indeed confirmed by the results for the stimulated 
emission signal reported below. 



C. Calculating the emission signal 

Next we calculate the time-dependent total stimulated 
emission signal. After the pump pulse at t = 0, the sys- 
tem is assumed to be in the excited state surface accord- 
ing to a properly chosen initial preparation. The probe 
pulse is then assumed to be much faster than typical sol- 
vent time scales such that the environmental influence 
can be neglected during the emission process itself. The 
time-dependent emission signal can thus be expressed in 
terms of the reduced density matrix directly before and 
after the application of the probe pulse. For a probe 
pulse centered at time t with duration 28, the energy 
E{t) emitted during the transition is 



E(i) = (H ) p ( t+S ) 
= TY{H [p(t- 



■ (Ho) P (t-S) 
5)- P (t-5)}} 



(4.4) 



G 



with the reduced density matrix p(t). Herein the influ- 
ence of the bath during the emission process has been ne- 
glected. Adopting a matrix representation for pit") with 
respect to the electronic states \G) and \E), we notice 
that p(t' <t-5) = \E)p E (t')(E\, since for t' < t-S, the 
wave packet is located on the excited state surface. For 
t' > t — S, however, V(t') causes a population of other 
matrix elements as well. Since we arc interested in the 



emission signal, the trace in Eq. (4.4) allows us to fo- 
cus only on the diagonal elements. Using second-order 
perturbation theory in V(t), they read 

p G (t + S) = Ui, GE (t + S,t-8)p E (t- 8) (4.5) 

x U^ EG (t + S,t-S) , 
p E {t + S) = Uo(t + 8,t-S)p E (t-8)U Q 1 (t + S,t-S) 

+ \u 2 {t + 6,t-S)p E (t- 8) Uo\t + S,t-6) 

+ H.c. 

Here Uk(t,t') denotes the appropriate matrix element of 
the fcth term of the Dyson expansion for the time evolu- 
tion operator under H(t), 



U (t,t') 

U h E G (t,t') 



£W) = --2 



e -%H E (t-t') 



dt 1 €*(t 1 )e-i HE( - t - t ^e-^ HG< - tl ~ t '^ , 



h 2 j t 



dh I dt 2 e(t 1 )e*(t 2 )e-K 



±H E (t-t 2 ) 



X e ~-k H G{t2-tl) e -^H E {tl-t') 



with Ui,GE(t,t') = —U\ EG (t,t'). After some algebra, we 
obtain the time-resolved total emission signal in the form 

oo / 

E(t) = - £ (n\r) (s\n)[2[ui (s + 1/2) +ui A ] (4.6) 

n,r,s— \ 

x Re{pf s (t - S)e l ^ r - S)it - S) 

x / dt'e(t') / dt"e*(t") 
Jt-s Jt-s 

x exp{i [uj G (n + 1/2) - uj (r + 1/2) - uj A ]t'} 
x exp{i [uj (s + 1/2) + lu a - u G {n + l/2)]t"}\ 
-u G {n + l/2)p E (t - 5 ) e ^o(r- s )(t-5) 



r t+5 

x I dt' e(t' 



1+6 



dt"e*(t") 

t-S Jt-8 

x exp{i [uj G (n + 1/2) - uj (r + 1/2) - uj A }t'} 
x exp{i [oj (s + 1/2) + uj a - oj G {n + l/2)]t"} 



where \n) and |r, s) denote the vibronic eigenstates of H G 
and He, respectively. 



In principle, the above analysis can straightforwardly 
be extended in order to incorporate the pump pulse. The 
resulting initial reduced density matrix is then composed 



of four different contributions, namely those in Eq. (4.5) 
and the two nondiagonal terms. The subsequent time 
evolution with both electronic surfaces coupled to the 
bath could then be treated in a similar way as presented 
in Sec. 0. 



D. Results for the reaction center 



Figure || shows the time-resolved stimulated emission 
signal for different probe wavelengths X pr at T = 10 K. At 
such a low temperature, the difference between the fac- 
torized and correlated preparation is very small and can 
hardly be resolved in Fig. ||. The qualitative features of 
the emission signal can be understood within the wave- 
packet picture by relating X pr to a particular value of 
the nuclear coordinate q, as is seen by plotting the dis- 
crete Fourier-transformed emission spectrum A(X pr ) at 
the frequency u> corresponding to the ground oscillation, 
see Fig. |^. This frequency, determined from the imag- 
inary part of the roots of Eq. (A12), is 97.2 cm -1 for 



SP II but deviates less than 1% from u>o for SP I. The 
maxima in A(X pr ) then correspond to the left (q = 0) 
and right (q = 2qo) turning points of the undamped nu- 
clear wave-packet, while the minimum is related to the 
bottom of the potential surface (q = qo) in Fig. ^. Due to 
the finite pulse duration and the different Franck-Condon 
overlap factors for q > q and q < qo, the correspond- 
ing value of A pr differs from 921.5 nm, particularly at 
high temperatures. Since the turning points are passed 
once per period but the bottom is visited twice, the cor- 
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FIG. 5. Emission signal E(t) [in arbitrary units] for dif- 
ferent probe wavelengths \ pr at T = 10 K for the factorized 
preparation. The solid (dashed) curve is for SP I (SP II). For 
clarity, curves for subsequent values of A pr have been shifted 
vertically. 
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responding emission signals should be oscillatory with 
frequency Q and 2u>, respectively JT^ |. This behavior is 
indeed found in Fig. ||. Focusing on SP I, the emission sig- 
nal at X pr = 1000 nm, corresponding to the right turning 
point, exhibits a phase shift of 7r and a smaller ampli- 
tude compared to X pr = 870 nm. This can be explained 
by noting that the right turning point is reached half a 
period later than the left one, whence the most signifi- 
cant initial contribution is damped more strongly. Apart 
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FIG. 6. Emission spectrum A(X pr ) [in arbitrary units] at 
u> for different temperatures. In (a) [(b)], we have taken SP I 
with a factorized [correlated] preparation. In (c) [(d)], we 
have taken SP II with a factorized [correlated] preparation. 

from the decay of the oscillations in E(t), damping of 
the nuclear motion is also reflected in a finite amplitude 
A(Xp r ) > at the minimum. For wavelengths X pr away 
from the turning points or the bottom, there are two dif- 
ferent time intervals between subsequent passings of the 
transition region. This leads to the splitting of the max- 
ima in the short-time emission signal observed in Fig. |^. 

With increasing temperature, according to our argu- 
mentation above, the influence of the initial preparation 
should become more and more important. This is seen in 
the temperature dependence of the spectrum and of the 
emission signal shown in Figures |^ and |7|, respectively. 
We first focus on the effects seen in FigTT(| The mini- 
mum of A(Xp r ) shifts towards higher wavelengths with 
increasing temperature. This can be rationalized by not- 
ing that for ujq < wg, the energy gap between the nth 
vibronic eigenstates of the excited and ground state de- 
creases with n, and that high-order eigenstates become 
more important in the spectral decomposition of p(t) at 
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FIG. 7. Emission signal E(t) for A pr = 950 nm and several 
temperatures for (a) SP I and (b) SP II. The solid (dashed) 
curve is for the correlated (factorized) preparation. Curves 
for subsequent temperatures have been shifted vertically by 
the same amount. 

higher temperatures. In fact, additional calculations for 
log — wo/2 [not shown here] yield a similar shift towards 
smaller X pr . This shift is more pronounced for the cor- 
related initial preparation, but depends only weakly on 
the spectral density of the environmental modes. 

Figure 7 shows the temperature dependence of the 
emission signal at X pr = 950 nm, corresponding to 
qo < q < 2qo, with q approaching the bottom q = qo 
as the temperature is increased. Furthermore, Figure 8 
shows that er(t) under the factorized preparation experi- 
ences a phase shift of almost ir at T > 100 K compared 
to the correlated preparation. Such an unphysical phase 
shift arises as a relict of the unperturbed evolution of a 
harmonic oscillator. Furthermore, since <r(0) = <t9 in- 
creases with temperature, a negative initial slope results 
under the factorized preparation, see Eq. ( p. 11 ). No- 
tably, for the correlated preparation, the maxima in E(t) 
stay always close to the equilibrium values, in marked 
contrast to the factorized preparation but in accordance 
with the experimental data of Ref. |l4j]. A similar be- 
havior is seen in the variance shown in Fig. |^. For the 
correlated preparation, the maxima in a(t) occur every 
(k + l/2)th period, corresponding to the passing of the 
bottom q — qo, and they are very close to their equilib- 
rium value erg. Therefore, for the correlated preparation, 
the independent expectation values (q(t)} and u{t) are al- 
ways close to their equilibrium values when passing the 
bottom q = qo. On the other hand, for the factorized 
preparation, the phase shift in a{t) results in a bunching 
of the wave packet when passing the bottom. This causes 
even qualitatively different emission signals. We conclude 
that at high temperatures several unphysical effects are 
introduced by using a factorized initial preparation, and 
a wrong description of the emission spectrum may result. 
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Temperature dependence of a(t) [in units of 



5 E = h/myjul + -yuj D } for (a) SP I and (b) SP II. The solid 
(dashed) curve is for the correlated (factorized) preparation. 



V. CONCLUSIONS 

In this work, we have formulated a dissipative wave- 
packet approach towards a detailed theoretical descrip- 
tion of stimulated emission pump-probe experiments un- 
der condensed-phase conditions. Assuming harmonic 
surfaces for both the ground and the excited state, the 
Gaussian nature of the wave packet describing the coher- 
ent nuclear motion allows for an exact treatment even 
if strong damping by environmental modes is present. 
Modelling the environmental modes by a set of infinitely 
many effective harmonic oscillators with a suitably cho- 
sen spectral density, it is then possible to make detailed 
predictions for the stimulated emission signal and for the 
corresponding spectra. While the spectral density is in 
principle accessible in terms of MD simulations, we have 
studied two model spectral densities in this work. A par- 
ticular advantage of our approach is the possibility of 
treating different initial preparations of the wavepacket- 
plus-bath complex directly after the pump pulse (t = 0) . 
A more realistic calculation should also explicitcly study 
the pump pulse, which can in principle be done along 



the same lines. Under such a formalism, a correct choice 
for the initial preparation of the wavepacket-plus-bath 
complex before the pump pulse will be important and is 
expected to lead to similar effects. 

The recent experiments by Vos et al. |l^-|l6[] on the 
bacterial photosynthctic reaction center have been an- 
alyzed using this formalism. Due to our assumptions 
about the pump pulse, the possibly important impul- 
sive resonant Raman contribution was not taken into ac- 
count here. While some of the qualitative features of the 
coherent nuclear motion have been discussed before us- 
ing simpler arguments JI3], our approach can allow for 
a fully quantum-mechanical comparison of experimental 
data with theory. Even in the absence of detailed knowl- 
edge about the environmental spectral density conclu- 
sions of relevance to the interpretation of experimental 
results can be extracted from our analysis. In particular, 
we have shown that at high temperatures, the assump- 
tion of a factorized initial state leads to large differences 
from the theoretical predictions under a more realistic 
correlated initial state. 

Finally it should be stressed that the approach pre- 
sented here can be applied to other pump-probe spec- 
troscopy setups as well. In particular, if the excited state 
surface is weakly coupled to another surface, as happens, 
e.g., in the primary electron transfer step in the reaction 
center, transitions to this surface are expected to modify 
the emission signal. A theoretical description of such a 
situation can be given in terms of spin-boson type models 
[E0|| and will be elaborated elsewhere. 
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This appendix contains the expressions for the vari- 
ances appearing in the general reduced density matrix 
fl3.7| ). Depending on the initial preparation, we get dif- 
ferent results as described below. 

For a factorized preparation, the three independent ex- 
pectation values are given byEqs. (|3.^3.iq). The various 
quantities appearing therein read as follows. For a spe- 
cific spectral density J(lo), the bath correlation function 
is 

1 f°° 

L(t) = — / du J(uj){coth(ujhf3/2) cos(u>i) — i sin(wt)} , 
7T Jo 

and the Laplace-transformed damping kernel reads, 

7(2) = / dw' V , 2j _ 2 • ( A1 ) 

7rm J uj uj + z 

The functions K q (t) and K p (t) are defined as 

K q (t) = / dt' G{t') f dt" G(t")L'(t' - t") , (A2) 
Jo Jo 

K p {t)= ( dt'G(t') f dt" G{t")L'{t' -t") , (A3) 
Jo Jo 
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with L'(t) = KeL(t) and G(t) being the inverse Laplace 
transform of 



G(z)= [z 2 + zj(z)+u J 2 ] 



(A4) 



In order to ensure that the wave packet relaxes to q = qa 
at long times, the bath must fulfill the condition G(t — > 
oo) = 0. Otherwise it merely leads to a mass renormal- 
ization but not to truly dissipative behavior. 



To obtain the time-dependent variances (3.9) and 
( |3 .10 ) in practice, it is necessary to find a more conve- 
nient form of Eqs. jA2| ) and JA3| ). By following Ref. 
we obtain [v n = 2"Kn/%(3) 



^) = Sf 2 E{( 1 - 6 " 1 ( i/ ») / dt'G{t')e-^' 

x f dt'G(t')e Vnt ' -G 2 (t)+ f dt' G(t')e- Vnt '\ 
Jo Jo ' 

- { ^2 - uj 2 J* dt' G{t')\ J* dt' G(t')}] (A5) 



and 



K p (t) = ^ [2 J2{ ("n ~ G-\v n ) dt' G{t')e- V ^ 



n=0 

t ft 

\ / dt' G(t')e Unt ' - i/„ / dt' G(t')e~ Vnt ' 
'0 Jo 

-G 2 (t) + l}+w 2 G 2 (i) . 
For i — > 00, the equilibrium variances then follow as 
2 



(A6) 



a l3 = a E 6 (KD' 

n— — 00 



■/?) 



E t 1 -^ 



(A7) 
(A8) 



n— — 00 



For the correlated preparation discussed in Sec. |IV B| , 
while (q(t)} stays the same as under the factorized prepa- 
ration, the variances now read 



<r(t) = — [A G G 2 (t) + Q G G 2 (t)+B(t) (A9) 
m V 

+ 2{A G G(t) G(t) C+(t) - G 2 (t) C+(t)} 

+ ±km) , 

([Ap} 2 )(t) = 7wti(a g G 2 (<) + tt G G 2 (i) + 5(t) (A10) 
+ 2\A G G(t)C 1 (t) - G(i)C 2 (i)j + 

Herein the various quantities are given as follows |l3| 

00 

nG = np E g g(KI)(<4 + KWKD), 



1 00 

a g=« H G G (K\), 

n—~oo 

B ^ = m E G G (Wn\) ds du [g n (s)g n (u) 

n— — oo 

-/n(»)/n(«M-«)G(t-tt), 
00 /"*/"* 
S ( t ) = hQ E G G(Wn\) ds du [g n {s)g n {u) 

n— — oo 

-/*(»)/n(tt)]G(t-«)G(t-tt), 



where G G is given by Eq. (A4) with cl>o being replaced 
by w G . The functions g n and /„ are given by 



Sn(s) 



1 



do; J(w) 



2w 



nor ,/q w 2 + f 2 



dio J{oj) 



2v„ 



cos(ws) , 



sin (us) 



Furthermore, we have used the abbreviations 



1 '.-V.' 

Cl{s)= hdT E G G (W n \)g n (s) 

n— — oo 
00 



C2{s) = ip E £ G (klK/„( s ) , 

G(< - s) 



n— — 00 



#(*) 
G(t) 



o dsCi{s) - G(t) ' 
t 

dsC<(«)G(t-a) . 



Next we briefly discuss the case of an ohmic bath with 
a Drude cutoff, see Eq. (3.4). The damping kernel then 
exhibits exponential decay, j(t) = jujd exp[— u>nt], with 
the Laplace transform 



7(z) 



COD + Z 



(All) 



Defining Aj for i — 1, 2, 3 as the roots of the cubic equa- 
tion 

z a — ujdz 2 + {"fujD + w 2 )z — luqLUd — , (A12) 
one obtains 

3 



G(z) 



Z + OJD 



E 



A; 



(z + Ai)(z + A 2 )(z + A 3 ) + \ 

where the coefficients A, follow as 
Xi(u>D - Aj) 



, (A13) 



Ai = 



2A?-u^(A?-u, 2 ) 



From Eq. ( A13[ ) we arrive at the simple result G(t) — 
J7 Aj exp[— Ait] . All variances can then be evaluated in 
closed form p2[. As the resulting expressions are very 
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lengthy but can be straightforwardly obtained by follow- 
ing the above steps, we refrain from quoting them here. 
To locate the coherent-to-incoherent transition, we 



note that the cubic equation (A12) has either three real 
solutions, or one real and two complex conjugate ones. 
In the latter case, G[t) and therefore (q(t)) will exhibit 
coherent oscillations. The critical value 7 C then follows 
from the condition D(fj,u>D) = 0, where 7 = 7/cJO) 
tin = wo/wo, and 

D ^u )D )=^ + f(±-^.) (A14) 
\uj d 4 J 

( 3 \ 1 2 

For arbitrary loq and u>d, there is exactly one positive 
value 7 = 7c solving the condition D = 0. 
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